
#######################################################################################

#load packages
library (Matrix); library (lme4); library (performance); library (lmerTest)
library (sjstats); library (stargazer); library (tidyverse); library("Hmisc")
library(car)

#get rid of scientific notation
options(scipen = 999)

#import data 

setwd ()

data = haven::read_dta("Processed data/merged_data_for_analysis-centered.dta") 

#KEEP ONLY KEY VARIABLES
data = data[ , c("neighborhood",
                 "neighborhood_key",
                 "week",
                 "nh_week",
                 "year",
                 "period1",
                 "period2",
                 "period3",
                 "vstop_wgt_avg_3wk_dev",
                 "pstop_wgt_avg_3wk_dev",
                 "arrest_drug_wgt_avg_3wk_dev",
                 "arrest_disorder_wgt_avg_3wk_dev",
                 "vstop_equal_3wk_dev",
                 "pstop_equal_3wk_dev",
                 "arrest_drug_equal_3wk_dev",
                 "arrest_disorder_equal_3wk_dev",
                 "vstop_3wk_dev19",
                 "pstop_3wk_dev19",
                 "arrest_drug_3wk_dev19",
                 "arrest_disorder_3wk_dev19",
                 "violent",
                 "crime_assault",
                 "crime_mvt",
                 "property",
                 "pstop",
                 "vstop",
                 "arrest_drug",
                 "arrest_disorder",
                 "accident",
                 "pct_black",
                 "pct_hispanic_any",
                 "prcp",
                 "tavg",
                 "AQI",
                 "total_pop",
                 "disadvantage",
                 "immigration",
                 "open_table")]

##eliminate observations pre-2020
data=data[data$year %in% 2020, ]
table(data$year)

#triple check that everything is centered appropriately
pastecs::stat.desc(data, basic=F)

###KEY VARIABLES###########################################################################
#mediators: vstop_wgt_avg_3wk_dev + arrest_drug_wgt_avg_3wk_dev + pstop_wgt_avg_3wk_dev + arrest_disorder_wgt_avg_3wk_dev
#controls w/ census factors-no dev: accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table
#between nhood effects: vstop_M pstop_M arrest_drug_M arrest_disorder_M
#within nhood effects: dev_vstop dev_pstop dev_arrest_drug dev_arrest_disorder 
######################################################################################################################################################################################
 


#check for spatial autocorrelation
library(sf); library(spdep)

s = st_read("Raw data/Neighborhood shapefiles/neighborhoods.shp")

#match shapefile names to merge file
s$neighborhood_key = str_replace(s$NBHD_NAME, " - ", "-")

data2 = dplyr::select(data,
                      neighborhood, neighborhood_key, 
                      week, nh_week, period2, period3, 
                      violent, property, crime_assault, crime_mvt)
shape = merge(s, data2, by = "neighborhood_key")
shape$NBHD_ID = shape$neighborhood #correct sorting issue in neighborhood numbers to match corrected merge file (neighborhood from Stata file, NBH_ID from shape)

#sort all datasets by nhweek (the stata nhweek) to avoid potential mismatching
data = data[order(data$neighborhood, data$nh_week), ]
shape = shape[order(shape$neighborhood, shape$nh_week), ]

#define neighboring polygons - going w/ queen
nb = poly2nb(shape, queen = T)

#get weights
lw = nb2listw(nb, style="W", zero.policy = T)

#Moran's I
moran_v = moran(shape$violent, lw, length(nb), Szero(lw))[1]
moran_v

moran_p = moran(shape$property, lw, length(nb), Szero(lw))[1]
moran_p

moran_a = moran(shape$crime_assault, lw, length(nb), Szero(lw))[1]
moran_a

moran_m = moran(shape$crime_mvt, lw, length(nb), Szero(lw))[1]
moran_m

#Test Moran's I
moran_v_test = moran.test(shape$violent, lw, alternative = "greater")
moran_v_test

moran_p_test = moran.test(shape$property, lw, alternative = "greater")
moran_p_test

moran_a_test = moran.test(shape$crime_assault, lw, alternative = "greater")
moran_a_test

moran_m_test = moran.test(shape$crime_mvt, lw, alternative = "greater")
moran_m_test

#create lagged measures by nhood - https://rpubs.com/erikaaldisa/spatialweights
violent_lag = lag.listw(lw, shape$violent)
lag.list.v = list(shape$NBHD_ID, shape$nh_week, lag.listw(lw, shape$violent))
lag.res.v = as.data.frame(lag.list.v)
colnames(lag.res.v) = c("NBHD_ID", "nh_week", "lag violent (queen)")
head(lag.res.v)
data = merge(data, lag.res.v, by = "nh_week")

property_lag = lag.listw(lw, shape$property)
lag.list.p = list(shape$NBHD_ID, shape$nh_week, lag.listw(lw, shape$property))
lag.res.p = as.data.frame(lag.list.p)
colnames(lag.res.p) = c("NBHD_ID", "nh_week", "lag property (queen)")
head(lag.res.p)
data = merge(data, lag.res.p, by = "nh_week")

assault_lag = lag.listw(lw, shape$crime_assault)
lag.list.a = list(shape$NBHD_ID, shape$nh_week, lag.listw(lw, shape$crime_assault))
lag.res.a = as.data.frame(lag.list.a)
colnames(lag.res.a) = c("NBHD_ID", "nh_week", "lag assault (queen)")
head(lag.res.a)
data = merge(data, lag.res.a, by = "nh_week")

mvt_lag = lag.listw(lw, shape$crime_mvt)
lag.list.m = list(shape$NBHD_ID, shape$nh_week, lag.listw(lw, shape$crime_mvt))
lag.res.m = as.data.frame(lag.list.m)
colnames(lag.res.m) = c("NBHD_ID", "nh_week", "lag mvt (queen)")
head(lag.res.m)
data = merge(data, lag.res.m, by = "nh_week")
#note: warning about duplicated NBHD_ID var from appending multiple - use to verify proper data order

save(data, file="data_w_lags.RData")




#PANEL 1 - OLS - MEDIATOR DEVIATIONS AS OUTCOMES###########################################################################################
#Pedestrian stops
t1p1_pstop = lmer (pstop_wgt_avg_3wk_dev ~ period2 + period3 + 
                  accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + 
                  (1 + period2 + period3|neighborhood), 
                data=data)
class(t1p1_pstop) = "lmerMod"
summary(t1p1_pstop)
BIC(t1p1_pstop)
AIC(t1p1_pstop)
logLik(t1p1_pstop, REML = F)

#Vehicle stops
t1p1_vstop = lmer (vstop_wgt_avg_3wk_dev ~ period2 + period3 + 
                  accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + 
                  (1 + period2 + period3|neighborhood), 
                data=data)
class(t1p1_vstop) = "lmerMod"
summary(t1p1_vstop)
BIC(t1p1_vstop)
AIC(t1p1_vstop)
logLik(t1p1_vstop, REML = F)

#Drug arrests
t1p1_darrest = lmer (arrest_drug_wgt_avg_3wk_dev ~ period2 + period3 + 
                    accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + 
                    (1 + period2 + period3|neighborhood), 
                  data=data)
class(t1p1_darrest) = "lmerMod"
summary(t1p1_darrest)  
BIC(t1p1_darrest)
AIC(t1p1_darrest)
logLik(t1p1_darrest, REML = F)

#Disorder arrests
t1p1_disarrest = lmer (arrest_disorder_wgt_avg_3wk_dev ~ period2 + period3 + 
                      accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + 
                      (1 + period2 + period3|neighborhood), 
                    data=data)
class(t1p1_disarrest) = "lmerMod"
summary(t1p1_disarrest) 
BIC(t1p1_disarrest)
AIC(t1p1_disarrest)
logLik(t1p1_disarrest, REML = F)

stargazer(t1p1_pstop, t1p1_vstop, t1p1_darrest, t1p1_disarrest, 
          type="html", star.cutoffs = c(0.05, 0.01, 0.001), out = "T1P1-Mediator deviations as outcomes.doc")





#PANEL 2 - NEGATIVE BINOMIAL - VIOLENCE COUNTS DIRECT EFFECTS
t1p2_violent = glmer.nb (violent ~ period2 + period3 + 
                                        accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag violent (queen)` + 
                                        (1 + period2 + period3|neighborhood), 
                                      data=data, nAGQ = 0)
summary(t1p2_violent)
BIC(t1p2_violent)
AIC(t1p2_violent)
logLik(t1p2_violent, REML = F)


###PERIODS + CONFOUNDERS + PSTOPS W/ VCS
t1p2_pstop = glmer.nb (violent ~ period2 + period3 + 
                                              pstop_wgt_avg_3wk_dev +
                                              accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag violent (queen)` +
                                              (1 + period2 + period3 + pstop_wgt_avg_3wk_dev|neighborhood), 
                                            data=data, family=poisson, nAGQ = 0)
summary(t1p2_pstop)
BIC(t1p2_pstop)
AIC(t1p2_pstop)
logLik(t1p2_pstop, REML = F)


###PERIODS + CONFOUNDERS + VSTOPS W/ VCS
t1p2_vstop = glmer.nb (violent ~ period2 + period3 + 
                                              vstop_wgt_avg_3wk_dev +
                                              accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag violent (queen)` +
                                              (1 + period2 + period3 + vstop_wgt_avg_3wk_dev|neighborhood), 
                                            data=data, family=poisson, nAGQ = 0)
summary(t1p2_vstop)
BIC(t1p2_vstop)
AIC(t1p2_vstop)
logLik(t1p2_vstop, REML = F)


##PERIODS + CONFOUNDERS + DARRESTS W/ VCS
t1p2_darrest = glmer.nb (violent ~ period2 + period3 + 
                                                arrest_drug_wgt_avg_3wk_dev +
                                                accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag violent (queen)` +
                                                (1 + period2 + period3 + arrest_drug_wgt_avg_3wk_dev|neighborhood), 
                                              data=data, family=poisson, nAGQ = 0)
summary(t1p2_darrest)
BIC(t1p2_darrest)
AIC(t1p2_darrest)
logLik(t1p2_darrest, REML = F)

###PERIODS + CONFOUNDERS + DISARRESTS W/ VCS
t1p2_disarrest = glmer.nb (violent ~ period2 + period3 + 
                                                  arrest_disorder_wgt_avg_3wk_dev +
                                                  accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag violent (queen)` +
                                                  (1 + period2 + period3 + arrest_disorder_wgt_avg_3wk_dev|neighborhood), 
                                                data=data, family=poisson, nAGQ = 0)
summary(t1p2_disarrest)
BIC(t1p2_disarrest)
AIC(t1p2_disarrest)
logLik(t1p2_disarrest, REML = F)


#get tables of fixed effects and add REs by hand
stargazer(t1p2_violent, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "T1P2-Violence counts periods and confounders - fixed only.doc")
stargazer(t1p2_pstop, t1p2_vstop, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "T1P2-Violence counts periods and confounders pstops and vstops - fixed only.doc")
stargazer(t1p2_darrest, t1p2_disarrest, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "T1P2-Violence counts periods and confounders drug and disorder arrests - fixed only.doc")

#PANEL 3 - PROPERTY COUNT MODELS - NEGATIVE BINOMIAL DIRECT EFFECTS###########################################################################################

###PERIODS + CONFROUNDERS W/ VCS

t1p3_property = glmer.nb (property ~ period2 + period3 + 
                                         accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag property (queen)` +
                                         (1 + period2 + period3|neighborhood), 
                                       data=data, nAGQ = 0)
summary(t1p3_property)
BIC(t1p3_property)
AIC(t1p3_property)
logLik(t1p3_property, REML = F)


###PERIODS + CONFOUNDERS + PSTOPS W/ VCS
t1p3_pstop = glmer.nb (property ~ period2 + period3 + 
                                               pstop_wgt_avg_3wk_dev +
                                               accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag property (queen)` + 
                                               (1 + period2 + period3 + pstop_wgt_avg_3wk_dev|neighborhood), 
                                             data=data, family=poisson, nAGQ = 0)
summary(t1p3_pstop)
BIC(t1p3_pstop)
AIC(t1p3_pstop)
logLik(t1p3_pstop, REML = F)



###PERIODS + CONFOUNDERS + VSTOPS W/ VCS
t1p3_vstop = glmer.nb (property ~ period2 + period3 + 
                                               vstop_wgt_avg_3wk_dev +
                                               accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag property (queen)` + 
                                               (1 + period2 + period3 + vstop_wgt_avg_3wk_dev|neighborhood), 
                                             data=data, family=poisson, nAGQ = 0)
summary(t1p3_vstop)
BIC(t1p3_vstop)
AIC(t1p3_vstop)
logLik(t1p3_vstop, REML = F)


###PERIODS + CONFOUNDERS + DARRESTS W/ VCS
t1p3_darrest = glmer.nb (property ~ period2 + period3 + 
                                                 arrest_drug_wgt_avg_3wk_dev +
                                                 accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag property (queen)` +
                                                 (1 + period2 + period3 + arrest_drug_wgt_avg_3wk_dev|neighborhood), 
                                               data=data, family=poisson, nAGQ = 0)
summary(t1p3_darrest)
BIC(t1p3_darrest)
AIC(t1p3_darrest)
logLik(t1p3_darrest, REML = F)


###PERIODS + CONFOUNDERS + DISARRESTS W/ VCS
t1p3_disarrest = glmer.nb (property ~ period2 + period3 + 
                                                   arrest_disorder_wgt_avg_3wk_dev +
                                                   accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag property (queen)` + 
                                                   (1 + period2 + period3 + arrest_disorder_wgt_avg_3wk_dev|neighborhood), 
                                                 data=data, family=poisson, nAGQ = 0)
summary(t1p3_disarrest)
BIC(t1p3_disarrest)
AIC(t1p3_disarrest)
logLik(t1p3_disarrest, REML = F)


#get tables of fixed effects and add REs by hand
stargazer(t1p3_property, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "T1P3-Property counts periods and confounders - fixed only.doc")
stargazer(t1p3_pstop, t1p3_vstop, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "T1P3-Property counts periods and confounders pstops and vstops - fixed only.doc")
stargazer(t1p3_darrest, t1p3_disarrest, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "T1P3-Property counts periods and confounders drug and disorder arrests - fixed only.doc")

save.image("denver_note-results-103024.RData")
#END OF MAIN MODELS############################################################################################################################










#SUPPLEMENTAL TABLES

#S4 - NEIGHBORHOODS ONLY
#PANEL 1 - OLS - MEDIATOR DEVIATIONS AS OUTCOMES###########################################################################################
#Pedestrian stops
st4p1_pstop = lmer (pstop_wgt_avg_3wk_dev ~ (1 |neighborhood), 
                   data=data)
class(st4p1_pstop) = "lmerMod"
summary(st4p1_pstop)
logLik(st4p1_pstop, REML = F)

#Vehicle stops
st4p1_vstop = lmer (vstop_wgt_avg_3wk_dev ~ (1 |neighborhood), 
                    data=data)
class(st4p1_vstop) = "lmerMod"
summary(st4p1_vstop)
logLik(st4p1_vstop, REML = F)

#Drug arrests
st4p1_darrest = lmer (arrest_drug_wgt_avg_3wk_dev ~ (1 |neighborhood), 
                    data=data)
class(st4p1_darrest) = "lmerMod"
summary(st4p1_darrest)
logLik(st4p1_darrest, REML = F)

#Disorder arrests
st4p1_disarrest = lmer (arrest_disorder_wgt_avg_3wk_dev ~ (1 |neighborhood), 
                      data=data)
class(st4p1_disarrest) = "lmerMod"
summary(st4p1_disarrest)
logLik(st4p1_disarrest, REML = F)

#PANEL 2 - NEGATIVE BINOMIAL - VIOLENCE AND PROPERTY COUNTS DIRECT EFFECTS
st4p2_violent = glmer.nb (violent ~ (1 |neighborhood), 
                         data=data, nAGQ = 0)
class(st4p2_violent) = "lmerMod"
summary(st4p2_violent)
logLik(st4p2_violent, REML = F)

st4p2_property = glmer.nb (property ~ (1 |neighborhood), 
                          data=data, nAGQ = 0)
class(st4p2_property) = "lmerMod"
summary(st4p2_property)
logLik(st4p2_property, REML = F)





#S8 - MVT MODELS - ###########################################################################################

###PERIODS + CONFROUNDERS W/ VCS
TS8P1_mvt = glmer.nb (crime_mvt ~ period2 + period3 + 
                                          accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag mvt (queen)` +
                                          (1 + period2 + period3|neighborhood), 
                                        data=data, nAGQ = 0)
summary(TS8P1_mvt)

###PERIODS + CONFOUNDERS + PSTOPS W/ VCS
TS8P1_mvt_pstop = glmer.nb (crime_mvt ~ period2 + period3 + 
                                                pstop_wgt_avg_3wk_dev +
                                                accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag mvt (queen)` +
                                                (1 + period2 + period3 + pstop_wgt_avg_3wk_dev|neighborhood), 
                                              data=data, nAGQ = 0)
summary(TS8P1_mvt_pstop)

###PERIODS + CONFOUNDERS + VSTOPS W/ VCS
TS8P1_mvt_vstop = glmer.nb (crime_mvt ~ period2 + period3 + 
                                                vstop_wgt_avg_3wk_dev +
                                                accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag mvt (queen)` +
                                                (1 + period2 + period3 + vstop_wgt_avg_3wk_dev|neighborhood), 
                                              data=data, nAGQ = 0)
summary(TS8P1_mvt_vstop)

###PERIODS + CONFOUNDERS + DARRESTS W/ VCS
TS8P1_mvt_darrest = glmer.nb (crime_mvt ~ period2 + period3 + 
                                                  arrest_drug_wgt_avg_3wk_dev +
                                                  accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag mvt (queen)` +
                                                  (1 + period2 + period3 + arrest_drug_wgt_avg_3wk_dev|neighborhood), 
                                                data=data, nAGQ = 0)
summary(TS8P1_mvt_darrest)

###PERIODS + CONFOUNDERS + DISARRESTS W/ VCS
TS8P1_mvt_disarrest = glmer.nb (crime_mvt ~ period2 + period3 + 
                                                    arrest_disorder_wgt_avg_3wk_dev +
                                                    accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag mvt (queen)` +
                                                    (1 + period2 + period3 + arrest_disorder_wgt_avg_3wk_dev|neighborhood), 
                                                  data=data, nAGQ = 0)
summary(TS8P1_mvt_disarrest)


#get tables of fixed effects and add REs by hand
stargazer(TS8P1_mvt, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "SuppMVT-Motor vehicle theft counts periods and confounders - fixed only.doc")
stargazer(TS8P1_mvt_pstop, TS8P1_mvt_vstop, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "SuppMVT-Motor vehicle theft counts periods and confounders pstops and vstops - fixed only.doc")
stargazer(TS8P1_mvt_darrest, TS8P1_mvt_disarrest, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "SuppMVT-Motor vehicle theft counts periods and confounders drug and disorder arrests - fixed only.doc")






#S9 - ASSAULT MODELS - ###########################################################################################
###PERIODS + CONFROUNDERS W/ VCS
TS9P1_assault = glmer.nb (crime_assault ~ period2 + period3 + 
                                              accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + assault_lag +
                                              (1 + period2 + period3|neighborhood), 
                                            data=data, nAGQ = 0)
summary(TS9P1_assault)

###PERIODS + CONFOUNDERS + PSTOPS W/ VCS
TS9P1_assault_pstop = glmer.nb (crime_assault ~ period2 + period3 + 
                                                    pstop_wgt_avg_3wk_dev +
                                                    accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + assault_lag +
                                                    (1 + period2 + period3 + pstop_wgt_avg_3wk_dev|neighborhood), 
                                                  data=data, nAGQ = 0)
summary(TS9P1_assault_pstop)

###PERIODS + CONFOUNDERS + VSTOPS W/ VCS
TS9P1_assault_vstop = glmer.nb (crime_assault ~ period2 + period3 + 
                                                    vstop_wgt_avg_3wk_dev +
                                                    accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + assault_lag +
                                                    (1 + period2 + period3 + vstop_wgt_avg_3wk_dev|neighborhood), 
                                                  data=data, nAGQ = 0)
summary(TS9P1_assault_vstop)

###PERIODS + CONFOUNDERS + DARRESTS W/ VCS
TS9P1_assault_darrest = glmer.nb (crime_assault ~ period2 + period3 + 
                                                      arrest_drug_wgt_avg_3wk_dev +
                                                      accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + assault_lag +
                                                      (1 + period2 + period3 + arrest_drug_wgt_avg_3wk_dev|neighborhood), 
                                                    data=data, nAGQ = 0)
summary(TS9P1_assault_darrest)

###PERIODS + CONFOUNDERS + DISARRESTS W/ VCS
TS9P1_assault_disarrest = glmer.nb (crime_assault ~ period2 + period3 + 
                                                         arrest_disorder_wgt_avg_3wk_dev +
                                                         accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + assault_lag +
                                                         (1 + period2 + period3 + arrest_disorder_wgt_avg_3wk_dev|neighborhood), 
                                                       data=data, nAGQ = 0)
summary(TS9P1_assault_disarrest)



#get tables of fixed effects and add REs by hand
stargazer(TS9P1_assault, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "SuppAssault-Assault counts periods and confounders - fixed only.doc")
stargazer(TS9P1_assault_pstop, TS9P1_assault_vstop, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "SuppAssault-Assault counts periods and confounders pstops and vstops - fixed only.doc")
stargazer(TS9P1_assault_darrest, TS9P1_assault_disarrest, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "SuppAssault-Assault counts periods and confounders drug and disorder arrests - fixed only.doc")







#SUPPLEMENTAL TABLES USING DIFFERENT DEVIATION WEIGHTS
#ST10 - EQUALLY WEIGHTED YEARS###########################################################################################################################################################
#PANEL 1 - OLS - MEDIATOR DEVIATIONS AS OUTCOMES###########################################################################################
#Pedestrian stops
TS10P1_pstop = lmer (pstop_equal_3wk_dev ~ period2 + period3 + 
                     accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + 
                     (1 + period2 + period3|neighborhood), 
                   data=data)
class(TS10P1_pstop) = "lmerMod"
summary(TS10P1_pstop)
BIC(TS10P1_pstop)

#Vehicle stops
TS10P1_vstop = lmer (vstop_equal_3wk_dev ~ period2 + period3 + 
                     accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + 
                     (1 + period2 + period3|neighborhood), 
                   data=data)
class(TS10P1_vstop) = "lmerMod"
summary(TS10P1_vstop)
BIC(TS10P1_vstop)

#Drug arrests
TS10P1_darrest = lmer (arrest_drug_equal_3wk_dev ~ period2 + period3 + 
                       accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + 
                       (1 + period2 + period3|neighborhood), 
                     data=data)
class(TS10P1_darrest) = "lmerMod"
summary(TS10P1_darrest)  
BIC(TS10P1_darrest)

#Disorder arrests
TS10P1_disarrest = lmer (arrest_disorder_equal_3wk_dev ~ period2 + period3 + 
                         accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + 
                         (1 + period2 + period3|neighborhood), 
                       data=data)
class(TS10P1_disarrest) = "lmerMod"
summary(TS10P1_disarrest) 
BIC(TS10P1_disarrest)


stargazer(TS10P1_pstop, TS10P1_vstop, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "ST2P1-Equal stop deviations as outcomes - fixed only.doc")

stargazer(TS10P1_darrest, TS10P1_disarrest, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "ST2P1-Equal arrest deviations as outcomes - fixed only.doc")

#PANEL 2 - NEGATIVE BINOMIAL - VIOLENCE COUNTS DIRECT EFFECTS
###PERIODS + CONFROUNDERS W/ VCS
TS10P2_violent = glmer.nb (violent ~ period2 + period3 + 
                           accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag violent (queen)` + 
                           (1 + period2 + period3|neighborhood), 
                         data=data, nAGQ = 0)
summary(TS10P2_violent)
class(TS10P2_violent) = "lmerMod"


###PERIODS + CONFOUNDERS + PSTOPS W/ VCS
TS10P2_violent_pstop = glmer.nb (violent ~ period2 + period3 + 
                                 pstop_equal_3wk_dev +
                                 accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag violent (queen)` +
                                 (1 + period2 + period3 + pstop_equal_3wk_dev|neighborhood), 
                               data=data, family=poisson, nAGQ = 0)
summary(TS10P2_violent_pstop)
class(TS10P2_violent_pstop) = "lmerMod"

###PERIODS + CONFOUNDERS + VSTOPS W/ VCS
TS10P2_violent_vstop = glmer.nb (violent ~ period2 + period3 + 
                                 vstop_equal_3wk_dev +
                                 accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag violent (queen)` +
                                 (1 + period2 + period3 + vstop_equal_3wk_dev|neighborhood), 
                               data=data, family=poisson, nAGQ = 0)
summary(TS10P2_violent_vstop)
class(TS10P2_violent_vstop) = "lmerMod"

###PERIODS + CONFOUNDERS + DARRESTS W/ VCS
TS10P2_violent_darrest = glmer.nb (violent ~ period2 + period3 + 
                                   arrest_drug_equal_3wk_dev +
                                   accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag violent (queen)` +
                                   (1 + period2 + period3 + arrest_drug_equal_3wk_dev|neighborhood), 
                                 data=data, family=poisson, nAGQ = 0)
summary(TS10P2_violent_darrest)
class(TS10P2_violent_darrest) = "lmerMod"


###PERIODS + CONFOUNDERS + DISARRESTS W/ VCS
TS10P2_violent_disarrest = glmer.nb (violent ~ period2 + period3 + 
                                     arrest_disorder_equal_3wk_dev +
                                     accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag violent (queen)` +
                                     (1 + period2 + period3 + arrest_disorder_equal_3wk_dev|neighborhood), 
                                   data=data, family=poisson, nAGQ = 0)
summary(TS10P2_violent_disarrest)
class(TS10P2_violent_disarrest) = "lmerMod"


#get tables of fixed effects and add REs by hand
stargazer(TS10P2_violent, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "ST2P2-Violence counts periods and confounders equal - fixed only.doc")
stargazer(TS10P2_violent_pstop, TS10P2_violent_vstop, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "ST2P2-Violence counts periods and confounders pstops and vstops equal - fixed only.doc")
stargazer(TS10P2_violent_darrest, TS10P2_violent_disarrest, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "ST2P2-Violence counts periods and confounders drug and disorder arrests equal - fixed only.doc")

#PANEL 3 - PROPERTY COUNT MODELS - NEGATIVE BINOMIAL DIRECT EFFECTS###########################################################################################
TS10P3_property = glmer.nb (property ~ period2 + period3 + 
                            accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag property (queen)` + 
                            (1 + period2 + period3 + pstop_equal_3wk_dev|neighborhood), 
                          data=data, family=poisson, nAGQ = 0)
summary(TS10P3_property)
class(TS10P3_property) = "lmerMod"


###PERIODS + CONFOUNDERS + PSTOPS W/ VCS
TS10P3_property_pstop = glmer.nb (property ~ period2 + period3 + 
                                  pstop_equal_3wk_dev +
                                  accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag property (queen)` + 
                                  (1 + period2 + period3 + pstop_equal_3wk_dev|neighborhood), 
                                data=data, family=poisson, nAGQ = 0)
summary(TS10P3_property_pstop)
class(TS10P3_property_pstop) = "lmerMod"


###PERIODS + CONFOUNDERS + VSTOPS W/ VCS
TS10P3_property_vstop = glmer.nb (property ~ period2 + period3 + 
                                  vstop_equal_3wk_dev +
                                  accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag property (queen)` + 
                                  (1 + period2 + period3 + vstop_equal_3wk_dev|neighborhood), 
                                data=data, family=poisson, nAGQ = 0)
summary(TS10P3_property_vstop)
class(TS10P3_property_vstop) = "lmerMod"



###PERIODS + CONFOUNDERS + DARRESTS W/ VCS
TS10P3_property_darrest = glmer.nb (property ~ period2 + period3 + 
                                    arrest_drug_equal_3wk_dev +
                                    accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag property (queen)` +
                                    (1 + period2 + period3 + arrest_drug_equal_3wk_dev|neighborhood), 
                                  data=data, family=poisson, nAGQ = 0)
summary(TS10P3_property_darrest)
class(TS10P3_property_darrest) = "lmerMod"


###PERIODS + CONFOUNDERS + DISARRESTS W/ VCS
TS10P3_property_disarrest = glmer.nb (property ~ period2 + period3 + 
                                      arrest_disorder_equal_3wk_dev +
                                      accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag property (queen)` + 
                                      (1 + period2 + period3 + arrest_disorder_equal_3wk_dev|neighborhood), 
                                    data=data, family=poisson, nAGQ = 0)
summary(TS10P3_property_disarrest)
class(TS10P3_property_disarrest) = "lmerMod"




#get tables of fixed effects and add REs by hand
stargazer(TS10P3_property, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "ST2P3-Property counts periods and confounders equal - fixed only.doc")
stargazer(TS10P3_property_pstop, TS10P3_property_vstop, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "ST2P3-Property counts periods and confounders pstops and vstops equal - fixed only.doc")
stargazer(TS10P3_property_darrest, TS10P3_property_disarrest, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "ST2P3-Property counts periods and confounders drug and disorder arrests equal - fixed only.doc")








#ALTERNATIVE DEVIATION SPECIFICATIONS
#S11 - 2019 ONLY########################################################################################################################################
#PANEL 1 - OLS - MEDIATOR DEVIATIONS AS OUTCOMES###########################################################################################
#Pedestrian stops
TS11P1_pstop = lmer (pstop_3wk_dev19 ~ period2 + period3 + 
                       accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + 
                       (1 + period2 + period3|neighborhood), 
                     data=data)
class(TS11P1_pstop) = "lmerMod"
summary(TS11P1_pstop)
BIC(TS11P1_pstop)

#Vehicle stops
TS11P1_vstop = lmer (vstop_3wk_dev19 ~ period2 + period3 + 
                       accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + 
                       (1 + period2 + period3|neighborhood), 
                     data=data)
class(TS11P1_vstop) = "lmerMod"
summary(TS11P1_vstop)
BIC(TS11P1_vstop)

#Drug arrests
TS11P1_darrest = lmer (arrest_drug_3wk_dev19 ~ period2 + period3 + 
                         accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + 
                         (1 + period2 + period3|neighborhood), 
                       data=data)
class(TS11P1_darrest) = "lmerMod"
summary(TS11P1_darrest)  
BIC(TS11P1_darrest)

#Disorder arrests
TS11P1_disarrest = lmer (arrest_disorder_3wk_dev19 ~ period2 + period3 + 
                           accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + 
                           (1 + period2 + period3|neighborhood), 
                         data=data)
class(TS11P1_disarrest) = "lmerMod"
summary(TS11P1_disarrest) 
BIC(TS11P1_disarrest)


stargazer(TS11P1_pstop, TS11P1_vstop, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "ST1P1-2019 Stop deviations as outcomes - fixed only.doc")

stargazer(TS11P1_darrest, TS11P1_disarrest, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "ST1P1-2019 Arrest deviations as outcomes - fixed only.doc")

#PANEL 2 - NEGATIVE BINOMIAL - VIOLENCE COUNTS DIRECT EFFECTS

###PERIODS + CONFROUNDERS W/ VCS
TS11P2_violent = glmer.nb (violent ~ period2 + period3 + 
                             accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag violent (queen)`  + 
                             (1 + period2 + period3|neighborhood), 
                           data=data, nAGQ = 0)
summary(TS11P2_violent)
class(TS11P2_violent) = "lmerMod"


###PERIODS + CONFOUNDERS + PSTOPS W/ VCS
TS11P2_violent_pstop = glmer.nb (violent ~ period2 + period3 + 
                                   pstop_3wk_dev19 +
                                   accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag violent (queen)`  +
                                   (1 + period2 + period3 + pstop_3wk_dev19|neighborhood), 
                                 data=data, family=poisson, nAGQ = 0)
summary(TS11P2_violent_pstop)
class(TS11P2_violent_pstop) = "lmerMod"

###PERIODS + CONFOUNDERS + VSTOPS W/ VCS
TS11P2_violent_vstop = glmer.nb (violent ~ period2 + period3 + 
                                   vstop_3wk_dev19 +
                                   accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag violent (queen)` +
                                   (1 + period2 + period3 + vstop_3wk_dev19|neighborhood), 
                                 data=data, family=poisson, nAGQ = 0)
summary(TS11P2_violent_vstop)
class(TS11P2_violent_vstop) = "lmerMod"

###PERIODS + CONFOUNDERS + DARRESTS W/ VCS
TS11P2_violent_darrest = glmer.nb (violent ~ period2 + period3 + 
                                     arrest_drug_3wk_dev19 +
                                     accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag violent (queen)` +
                                     (1 + period2 + period3 + arrest_drug_3wk_dev19|neighborhood), 
                                   data=data, family=poisson, nAGQ = 0)
summary(TS11P2_violent_darrest )
class(TS11P2_violent_darrest ) = "lmerMod"


###PERIODS + CONFOUNDERS + DISARRESTS W/ VCS
TS11P2_violent_disarrest = glmer.nb (violent ~ period2 + period3 + 
                                       arrest_disorder_3wk_dev19 +
                                       accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag violent (queen)` +
                                       (1 + period2 + period3 + arrest_disorder_3wk_dev19|neighborhood), 
                                     data=data, family=poisson, nAGQ = 0)
summary(TS11P2_violent_disarrest)
class(TS11P2_violent_disarrest) = "lmerMod"


#get tables of fixed effects and add REs by hand
stargazer(TS11P2_violent, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "ST1P2-Violence counts periods and confounders 2019 - fixed only.doc")
stargazer(TS11P2_violent_pstop, TS11P2_violent_vstop, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "ST1P2-Violence counts periods and confounders pstops and vstops 2019 - fixed only.doc")
stargazer(TS11P2_violent_darrest, TS11P2_violent_disarrest, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "ST1P2-Violence counts periods and confounders drug and disorder arrests 2019 - fixed only.doc")

#PANEL 3 - PROPERTY COUNT MODELS - NEGATIVE BINOMIAL DIRECT EFFECTS###########################################################################################

###PERIODS + CONFROUNDERS W/ VCS
TS11P3_property = glmer.nb (property ~ period2 + period3 + 
                              accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag property (queen)` +
                              (1 + period2 + period3|neighborhood), 
                            data=data, nAGQ = 0)
summary(TS11P3_property)
class(TS11P3_property) = "lmerMod"



###PERIODS + CONFOUNDERS + PSTOPS W/ VCS
TS11P3_property_pstop = glmer.nb (property ~ period2 + period3 + 
                                    pstop_3wk_dev19 +
                                    accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag property (queen)` + 
                                    (1 + period2 + period3 + pstop_3wk_dev19|neighborhood), 
                                  data=data, family=poisson, nAGQ = 0)
summary(TS11P3_property_pstop)
class(TS11P3_property_pstop) = "lmerMod"


###PERIODS + CONFOUNDERS + VSTOPS W/ VCS
TS11P3_property_vstop = glmer.nb (property ~ period2 + period3 + 
                                    vstop_3wk_dev19 +
                                    accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag property (queen)` +
                                    (1 + period2 + period3 + vstop_3wk_dev19|neighborhood), 
                                  data=data, family=poisson, nAGQ = 0)
summary(TS11P3_property_vstop)
class(TS11P3_property_vstop) = "lmerMod"



###PERIODS + CONFOUNDERS + DARRESTS W/ VCS
TS11P3_property_darrest = glmer.nb (property ~ period2 + period3 + 
                                      arrest_drug_3wk_dev19 +
                                      accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag property (queen)` +
                                      (1 + period2 + period3 + arrest_drug_3wk_dev19|neighborhood), 
                                    data=data, family=poisson, nAGQ = 0)
summary(TS11P3_property_darrest)
class(TS11P3_property_darrest) = "lmerMod"


###PERIODS + CONFOUNDERS + DISARRESTS W/ VCS
TS11P3_property_disarrest = glmer.nb (property ~ period2 + period3 + 
                                        arrest_disorder_3wk_dev19 +
                                        accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag property (queen)` +
                                        (1 + period2 + period3 + arrest_disorder_3wk_dev19|neighborhood), 
                                      data=data, family=poisson, nAGQ = 0)
summary(TS11P3_property_disarrest)
class(TS11P3_property_disarrest) = "lmerMod"




#get tables of fixed effects and add REs by hand
stargazer(TS11P3_property, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "ST1P3-Property counts periods and confounders 2019 - fixed only.doc")
stargazer(TS11P3_property_pstop, TS11P3_property_vstop, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "ST1P3-Property counts periods and confounders pstops and vstops 2019 - fixed only.doc")
stargazer(TS11P3_property_darrest, TS11P3_property_disarrest, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "ST1P3-Property counts periods and confounders drug and disorder arrests 2019 - fixed only.doc")






#s12 - TOTAL EFFECTS MODEL W/ COVID AND FLOYD PERIOD COMBINED##################################
data$period4=data$period1==0

#PANEL 1 - OLS - MEDIATOR DEVIATIONS AS OUTCOMES###########################################################################################
#Pedestrian stops
TS12P1_pstop = lmer (pstop_wgt_avg_3wk_dev ~ period4 + 
                            accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + 
                            (1 + period4|neighborhood), 
                          data=data)
class(TS12P1_pstop) = "lmerMod"
summary(TS12P1_pstop)
BIC(TS12P1_pstop)
AIC(TS12P1_pstop)
logLik(TS12P1_pstop, REML = F)

#Vehicle stops
TS12P1_vstop = lmer (vstop_wgt_avg_3wk_dev ~ period4 + 
                            accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + 
                            (1 + period4|neighborhood), 
                          data=data)
class(TS12P1_vstop) = "lmerMod"
summary(TS12P1_vstop)
BIC(TS12P1_vstop)
AIC(TS12P1_vstop)
logLik(TS12P1_vstop, REML = F)

#Drug arrests
TS12P1_darrest = lmer (arrest_drug_wgt_avg_3wk_dev ~ period4 + 
                              accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + 
                              (1 + period4|neighborhood), 
                            data=data)
class(TS12P1_darrest) = "lmerMod"
summary(TS12P1_darrest)  
BIC(TS12P1_darrest)
AIC(TS12P1_darrest)
logLik(TS12P1_darrest, REML = F)

#Disorder arrests
TS12P1_disarrest = lmer (arrest_disorder_wgt_avg_3wk_dev ~ period4 + 
                                accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + 
                                (1 + period4|neighborhood), 
                              data=data)
class(TS12P1_disarrest) = "lmerMod"
summary(TS12P1_disarrest) 
BIC(TS12P1_disarrest)
AIC(TS12P1_disarrest)
logLik(TS12P1_disarrest, REML = F)


#get tables of fixed effects and add REs by hand
stargazer(TS12P1_pstop, TS12P1_vstop, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "ST-Stop deviations total effect - fixed only.doc")
stargazer(TS12P1_darrest, TS12P1_disarrest, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "ST-Arrest deviations total effect - fixed only.doc")
###################################################################

#PANEL 2 - NEGATIVE BINOMIAL - VIOLENCE COUNTS DIRECT EFFECTS
TS12P2_violent = glmer.nb (violent ~ period4 + 
                                        accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + violent_lag + 
                                        (1 + period4|neighborhood), 
                                      data=data, nAGQ = 0)
summary(TS12P2_violent)


TS12P2_violent_pstop = glmer.nb (violent ~ period4 + 
                                              pstop_wgt_avg_3wk_dev +
                                              accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + violent_lag +
                                              (1 + period4 + pstop_wgt_avg_3wk_dev|neighborhood), 
                                            data=data, family=poisson, nAGQ = 0)
summary(TS12P2_violent_pstop)


TS12P2_violent_vstop = glmer.nb (violent ~ period4 + 
                                              vstop_wgt_avg_3wk_dev +
                                              accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + violent_lag +
                                              (1 + period4 + vstop_wgt_avg_3wk_dev|neighborhood), 
                                            data=data, family=poisson, nAGQ = 0)
summary(TS12P2_violent_vstop)


TS12P2_violent_darrest = glmer.nb (violent ~ period4 + 
                                                arrest_drug_wgt_avg_3wk_dev +
                                                accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + violent_lag +
                                                (1 + period4 + arrest_drug_wgt_avg_3wk_dev|neighborhood), 
                                              data=data, family=poisson, nAGQ = 0)
summary(TS12P2_violent_darrest)


TS12P2_violent_disarrest = glmer.nb (violent ~ period4 + 
                                                  arrest_disorder_wgt_avg_3wk_dev +
                                                  accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + violent_lag +
                                                  (1 + period4 + arrest_disorder_wgt_avg_3wk_dev|neighborhood), 
                                                data=data, family=poisson, nAGQ = 0)
summary(TS12P2_violent_disarrest)


#get tables of fixed effects and add REs by hand
stargazer(TS12P2_violent, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "ST-Violence counts periods and confounders total effect - fixed only.doc")
stargazer(TS12P2_violent_pstop, TS12P2_violent_vstop, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "ST-Violence counts periods and confounders pstops and vstops total effect - fixed only.doc")
stargazer(TS12P2_violent_darrest, TS12P2_violent_disarrest, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "ST-Violence counts periods and confounders drug and disorder arrests total effect - fixed only.doc")

#PANEL 3 - PROPERTY COUNT MODELS - NEGATIVE BINOMIAL DIRECT EFFECTS###########################################################################################

###PERIODS + CONFROUNDERS W/ VCS
TS12P3_property = glmer.nb (property ~ period4 + 
                                         accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + property_lag +
                                         (1 + period4|neighborhood), 
                                       data=data, nAGQ = 0)
summary(TS12P3_property)



TS12P3_property_pstop = glmer.nb (property ~ period4 + 
                                               pstop_wgt_avg_3wk_dev +
                                               accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + property_lag + 
                                               (1 + period4 + pstop_wgt_avg_3wk_dev|neighborhood), 
                                             data=data, family=poisson, nAGQ = 0)
summary(TS12P3_property_pstop)


TS12P3_property_vstop = glmer.nb (property ~ period4 + 
                                               vstop_wgt_avg_3wk_dev +
                                               accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + property_lag + 
                                               (1 + period4 + vstop_wgt_avg_3wk_dev|neighborhood), 
                                             data=data, family=poisson, nAGQ = 0)
summary(TS12P3_property_vstop)


TS12P3_property_darrest = glmer.nb (property ~ period4 + 
                                                 arrest_drug_wgt_avg_3wk_dev +
                                                 accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + property_lag +
                                                 (1 + period4 + arrest_drug_wgt_avg_3wk_dev|neighborhood), 
                                               data=data, family=poisson, nAGQ = 0)
summary(TS12P3_property_darrest)


TS12P3_property_disarrest = glmer.nb (property ~ period4 + 
                                                   arrest_disorder_wgt_avg_3wk_dev +
                                                   accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + property_lag + 
                                                   (1 + period4 + arrest_disorder_wgt_avg_3wk_dev|neighborhood), 
                                                 data=data, family=poisson, nAGQ = 0)
summary(TS12P3_property_disarrest)



#get tables of fixed effects and add REs by hand
stargazer(TS12P3_property, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "ST-Property counts periods and confounders total effect - fixed only.doc")
stargazer(TS12P3_property_pstop, TS12P3_property_vstop, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "ST-Property counts periods and confounders pstops and vstops total effect - fixed only.doc")
stargazer(TS12P3_property_darrest, TS12P3_property_disarrest, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "ST-Property counts periods and confounders drug and disorder arrests total effect - fixed only.doc")
################################################################################################################################################





#S13 - MODELS USING COVID ONLY AS PERIOD 0 AND FLOYD AS PERIOD 1#############################################################################################################
data$periodrr = NA
data$periodrr = ifelse(data$period2==1, 0, data$periodrr)
data$periodrr = ifelse(data$period3==1, 1, data$periodrr)

#PANEL 1 - OLS - MEDIATOR DEVIATIONS AS OUTCOMES###########################################################################################
#Pedestrian stops
TS13P1_pstop = lmer (pstop_wgt_avg_3wk_dev ~ periodrr + 
                            accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + 
                            (1 + periodrr|neighborhood), 
                          data=data)
class(TS13P1_pstop) = "lmerMod"
summary(TS13P1_pstop)
BIC(TS13P1_pstop)
AIC(TS13P1_pstop)
logLik(TS13P1_pstop, REML = F)

#Vehicle stops
TS13P1_vstop = lmer (vstop_wgt_avg_3wk_dev ~ periodrr + 
                            accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + 
                            (1 + periodrr|neighborhood), 
                          data=data)
class(TS13P1_vstop) = "lmerMod"
summary(TS13P1_vstop)
BIC(TS13P1_vstop)
AIC(TS13P1_vstop)
logLik(TS13P1_vstop, REML = F)

#Drug arrests
TS13P1_darrest = lmer (arrest_drug_wgt_avg_3wk_dev ~ periodrr + 
                              accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + 
                              (1 + periodrr|neighborhood), 
                            data=data)
class(TS13P1_darrest) = "lmerMod"
summary(TS13P1_darrest)  
BIC(TS13P1_darrest)
AIC(TS13P1_darrest)
logLik(TS13P1_darrest, REML = F)

#Disorder arrests
TS13P1_disarrest = lmer (arrest_disorder_wgt_avg_3wk_dev ~ periodrr + 
                                accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + 
                                (1 + periodrr|neighborhood), 
                              data=data)
class(TS13P1_disarrest) = "lmerMod"
summary(TS13P1_disarrest) 
BIC(TS13P1_disarrest)
AIC(TS13P1_disarrest)
logLik(TS13P1_disarrest, REML = F)


stargazer(TS13P1_pstop, TS13P1_vstop, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "T1P1-Stops as outcomes - fixed only.doc")

stargazer(TS13P1_darrest, TS13P1_disarrest, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "T1P1-Drug and disorder arrests as outcomes - fixed only.doc")




#PANEL 2 - NEGATIVE BINOMIAL - VIOLENCE COUNTS DIRECT EFFECTS
TS13P2_violent = glmer.nb (violent ~ periodrr + 
                                        accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag violent (queen)` + 
                                        (1 + periodrr|neighborhood), 
                                      data=data, nAGQ = 0)
summary(TS13P2_violent)


TS13P2_violent_pstop = glmer.nb (violent ~ periodrr + 
                                              pstop_wgt_avg_3wk_dev +
                                              accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag violent (queen)` +
                                              (1 + periodrr + pstop_wgt_avg_3wk_dev|neighborhood), 
                                            data=data, family=poisson, nAGQ = 0)
summary(TS13P2_violent_pstop)



TS13P2_violent_vstop = glmer.nb (violent ~ periodrr + 
                                              vstop_wgt_avg_3wk_dev +
                                              accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag violent (queen)` +
                                              (1 + periodrr + vstop_wgt_avg_3wk_dev|neighborhood), 
                                            data=data, family=poisson, nAGQ = 0)
summary(TS13P2_violent_vstop)


TS13P2_violent_darrest = glmer.nb (violent ~ periodrr + 
                                                arrest_drug_wgt_avg_3wk_dev +
                                                accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag violent (queen)` +
                                                (1 + periodrr + arrest_drug_wgt_avg_3wk_dev|neighborhood), 
                                              data=data, family=poisson, nAGQ = 0)
summary(TS13P2_violent_darrest)



TS13P2_violent_disarrest = glmer.nb (violent ~ periodrr + 
                                                  arrest_disorder_wgt_avg_3wk_dev +
                                                  accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag violent (queen)` +
                                                  (1 + periodrr + arrest_disorder_wgt_avg_3wk_dev|neighborhood), 
                                                data=data, family=poisson, nAGQ = 0)
summary(TS13P2_violent_disarrest)


#get tables of fixed effects and add REs by hand
stargazer(TS13P2_violent, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "T1P2-Violence counts periods and confounders - fixed only.doc")
stargazer(TS13P2_violent_pstop, TS13P2_violent_vstop, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "T1P2-Violence counts periods and confounders pstops and vstops - fixed only.doc")
stargazer(TS13P2_violent_darrest, TS13P2_violent_disarrest, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "T1P2-Violence counts periods and confounders drug and disorder arrests - fixed only.doc")

#PANEL 3 - PROPERTY COUNT MODELS - NEGATIVE BINOMIAL DIRECT EFFECTS###########################################################################################
TS13P3_property = glmer.nb (property ~ periodrr + 
                                         accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag property (queen)` +
                                         (1 + periodrr|neighborhood), 
                                       data=data, nAGQ = 0)
summary(TS13P3_property)

TS13P3_property_pstop = glmer.nb (property ~ periodrr + 
                                               pstop_wgt_avg_3wk_dev +
                                               accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag property (queen)` + 
                                               (1 + periodrr + pstop_wgt_avg_3wk_dev|neighborhood), 
                                             data=data, family=poisson, nAGQ = 0)
summary(TS13P3_property_pstop)



TS13P3_property_vstop = glmer.nb (property ~ periodrr + 
                                               vstop_wgt_avg_3wk_dev +
                                               accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag property (queen)` + 
                                               (1 + periodrr + vstop_wgt_avg_3wk_dev|neighborhood), 
                                             data=data, family=poisson, nAGQ = 0)
summary(TS13P3_property_vstop)



TS13P3_property_darrest = glmer.nb (property ~ periodrr + 
                                                 arrest_drug_wgt_avg_3wk_dev +
                                                 accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag property (queen)` +
                                                 (1 + periodrr + arrest_drug_wgt_avg_3wk_dev|neighborhood), 
                                               data=data, family=poisson, nAGQ = 0)
summary(TS13P3_property_darrest)



TS13P3_property_disarrest = glmer.nb (property ~ periodrr + 
                                                   arrest_disorder_wgt_avg_3wk_dev +
                                                   accident + total_pop + disadvantage + pct_black + pct_hispanic_any + immigration + prcp + tavg + AQI + open_table + `lag property (queen)` + 
                                                   (1 + periodrr + arrest_disorder_wgt_avg_3wk_dev|neighborhood), 
                                                 data=data, family=poisson, nAGQ = 0)
summary(TS13P3_property_disarrest)



#get tables of fixed effects and add REs by hand
stargazer(TS13P3_property, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "T1P3-Property counts periods and confounders - fixed only.doc")
stargazer(TS13P3_property_pstop, TS13P3_property_vstop, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "T1P3-Property counts periods and confounders pstops and vstops - fixed only.doc")
stargazer(TS13P3_property_darrest, TS13P3_property_disarrest, type="html", 
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "T1P3-Property counts periods and confounders drug and disorder arrests - fixed only.doc")

#END OF REVISED PERIOD MODELS############################################################################################################################
